A Model for the Elasticity of Compressed Emulsions 



On 



(N 



> 

O 
O 

m 
o 

OS 



X 



Martin-D. Lacasse 1 , Gary S. Grest 1 , Dov Levinc 2 , T.G. Mason 1 !] and D.A. Wcitz 1 ! 

Corporate Research Science Laboratories 
Exxon Research and Engineering Co., Annandale, NJ 08801 
2 Department of Physics, Technion, Haifa, 32000 Israel 
(December 13, 1995) 

We present a new model to describe the unusual elastic properties of compressed emulsions. The 
response of a single droplet under compression is investigated numerically for different Wigner-Seitz 
cells. The response is softer than harmonic, and depends on the coordination number of the droplet. 
Using these results, we propose a new effective inter-droplet potential which is used to determine the 
elastic response of a monodisperse collection of disordered droplets as a function of volume fraction. 
Our results are in excellent agreement with recent experiments. This suggests that anharmonicity, 
together with disorder, are responsible for the quasi-linear increase of G and fl observed at tp c . 



PACS numbers: 82.70.Kj, 81.40. Jj, 62.20.Dc 

Emulsions are materials with highly unusual elastic 
properties. They consist of droplets of one fluid dis- 
persed in a second fluid, with interfaces stabilized by a 
surfactant. Despite being comprised solely of fluids, they 
can be elastic solids, when droplets are compressed to 
a large volume fraction, tp, by an osmotic pressure, II. 
The origin of this elasticity is the intcrfacial energy of 
the droplets. At low volume fractions, their surface ten- 
sion, cr, ensures that the droplets are spherical in shape; 
however, at higher tp, the packing constraints force the 
droplet shapes to deform, storing energy. The applica- 
tion of a shear strain to a compressed emulsion causes 
the droplets to deform further, increasing their surface 
area, thereby storing elastic energy 




FIG. 1. Experimental values of the osmotic pressure (0) 
and elastic shear modulus (•) of different emulsions with 
monodisperse droplets. The data are normalized by o-/R, 
allowing data from emulsions with different droplet sizes to 
be compared. Results predicted by the present model for the 
shear modulus (+) and the osmotic pressure (line) are also 
shown. 



The experimentally measured |6j,|7| tp dependence of 
the static shear modulus G is remarkable in several ways. 
Typical data j?J obtained for monodisperse || silicone- 
oil-in-water emulsions are summarized in Fig. |l|. The 
solid symbols show G measured as a function of tp for 



several monodisperse emulsions of different droplet sizes. 
The data are normalized by a/R where R is the ra- 
dius of an undeformed droplet. The observed scaling 
by the Laplace pressure (2a /R) confirms the essential 
role of the interfacial energy. The magnitude of the 
scaled modulus increases approximately linearly, vary- 
ing as G ~ tp(tp — tp c ), where tp c is the volume fraction 
where the droplets are first deformed. For monodisperse 
droplets, tp c « 0.64, the maximum volume fraction that 
monodisperse spheres can be randomly packed for 
polydisperse droplets tp c is larger jfj), reflecting more ef- 
ficient packing. Surprisingly, the osmotic pressure re- 
quired to compress the emulsion is very similar to the 
shear modulus Jt],[io|], as shown by the open circles in 
Fig. [l]. This implies that the corresponding longitudi- 
nal bulk elastic modulus, or the bulk osmotic modulus, 
K = tpdll/dip, must differ significantly from the shear 
modulus as it must have a much sharper onset at tp c [Q. 

The behavior of the shear modulus of emulsions has 
been the subject of various theoretical studies. A first 
approach consists in considering a single droplet, and to 
obtain the properties of a periodic structure, assuming 
the strain remains affine [0J^J^]. While this can be done 
analytically in two dimensions, in three dimensions it can 
only be done through approximations. Results obtained 
from this approach consistently predict a sharp rise of G 
at tp c , followed by a much slower increase at larger tp. 
This suggests that the effects of disorder may be respon- 
sible for the quasi-linear rise of G at tp c found experimen- 
tally. The response to shear of disordered arrangements 
of droplets has also been investigated. However, due to 
the additional complexity this imposes, these calculations 
have been restricted to two dimensions [||,|],[ll]|lj] , where 
disorder needs to be introduced through polydispersity. 
Nevertheless, these two-dimensional models show a dif- 
ferent behavior than that observed experimentally for 
disordered, three-dimensional emulsions. In particular, 
two-dimensional systems do not exhibit a smooth nearly 
linear increase of G near tp c , although a smooth quadratic 
increase of II near tp c has recently been observed [jllj . 
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In this Letter, we present the first three-dimensional 
computer simulation modeling the elasticity of a disor- 
dered emulsion. In addition, we propose a new, more 
realistic inter-droplet potential. It is based on numer- 
ical results obtained by calculating the change in sur- 
face energy of a single droplet as it is compressed using 
Brakke's Surface Evolver (SE) software Our model 
captures the essential physics of the droplet interactions 
and their disordered, glassy packing. The droplets are 
soft, monodisperse spheres interacting through a purely 
repulsive potential. The droplet positions are allowed to 
relax upon application of shear. We calculate the osmotic 
pressure and the shear modulus as a function of volume 
fraction. We find that both the positional relaxation of 
the droplets and the form of the potential are essential 
in reproducing the experimental behavior of the shear 
modulus. 

We first investigate in more detail the interaction be- 
tween droplets. When two droplets are forced together, 
their spherical shapes are deformed and their surfaces 
develop flattened facets at contact. As a simplification, 
consider a droplet of radius R compressed between two 
parallel planes, each one located at a distance h from the 
center of the droplet. Naively, for small deformations 
where £ = (R — h)/ R is a dimensionless measure of com- 
pression, the resultant force, F, on the flat facets can be 
estimated by assuming that the radius of the droplet, and 
hence the Laplace pressure, remains unchanged. We thus 
have F w (2<j/R)8S, where SS is the area of the flattened 
facet. To linear order in the deformation, SS « 2irR 2 S£,, 
so that F » AircrRS^. Thus, the interaction between the 
facets of two droplets is often taken as a strictly repulsive 
harmonic spring of spring constant Ana acting between 
the centers of the spheres ||,|l2| . 

The harmonic-spring potential, while appealing, ig- 
nores the details of the response of the shape to defor- 
mation 0,0, and the possibility of coupling between 
the different facets on each droplet. Thus, to determine 
an improved potential, we investigated numerically the 
shape of a single droplet confined within space-filling 
polyhedral cells. Using SE we calculated the excess sur- 
face as the confinement is increased, under the constraint 
of a fixed droplet volume. As confining cells we investi- 
gated a rhombic dodecahedron (f.c.c), a truncated octa- 
hedron (b.c.c.) and a simple cube (s.c). In all cases, the 
compressions leave the center of mass unchanged, so that 
a distance h from the center can be defined. 

For simplicity, we shall define a — 1, and R = 1. For 
a droplet of total surface area A, we can thus define the 
excess energy as E\ = (A—4ir). Figure |^ shows the calcu- 
lated dependence of the excess energy per facet, E\/n, on 
£ for different types of cells, where n represents the num- 
ber of facets. The data are well described by [(R/h) 3 — l] a 
as shown by the dashed lines. For small £, this form re- 
duces to a simple power law, (3£) a , as is evidenced in 
Fig.H by the asymptotically linear behavior. This func- 
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FIG. 2. Excess surface energy excess per facet under 
different Wigner-Seitz constraining cells. The lines are 
obtained from a fit to Ei/n = C(<p/ip c — l) a , where 
ip/ifc = (R/h) 3 — 1/(1 — £) 3 . Curves are, from top to bottom: 
f.c.c. (n = 12, tp c = 7i"v2/6), b.c.c. (n = 8, ip c = 7r\/3/8), and 
s.c. (n — 6, ip c — 7r/6). Values of C and a are accurate to 
0.05. 

tional form is also equivalent to (ip — <p c ) a for space- filling 
structures. Fitting to Ei/n = C[(R/h) 3 - 1] Q , we find 
that both the coefficient, C, and the exponent, a, depend 
on n, as is evident from the data in Fig|^. Their values are 
shown in the table inserted in Fig. ^. The data selected 
for the fits lie in the £ interval 0-7%, the upper value 
corresponding to <p = 0.92 for an f.c.c. lattice. For the 
b.c.c. lattice we use n = 8 since second neighbors do not 
contribute to the energy for tp < 0.90 [H. Our results 
show explicitly that the response of a droplet to compres- 
sion is a non-local phenomenon: the response depends on 
the number of planes used to compress the droplet. One 
might expect that it also depends on the relative ori- 
entation of these planes. However, comparison between 
results obtained for iso-n configurations, for instance, the 
pentagon dodecahedron and the rhombic dodecahedron, 
suggests that the distribution of the compressing planes 
on the droplet surface has only a minor effect 0. Thus, 
we expect a similar functional form for E\, even when 
the planes are more randomly distributed as they would 
in a disordered emulsion. 

This behavior is different from the one encountered 
in two-dimensional systems. A minimum free surface is 
characterized by a uniform pressure or, equivalently by a 
uniform mean curvature. In two dimensions, the surface 
is parameterized by only one radius of curvature and the 
minimum free surface is always an arc of a circle. In 
this case, a harmonic potential is a good approximation. 
By contrast, in three dimensions, the response is always 
softer than harmonic since a(n) > 2 for all values of 
n investigated here; thus, the effective spring constant 
goes to zero as the distortion vanishes. We note that, as 
proposed, a n-dependent power-law leads to unphysical 
behavior for very small compressions, as the energies for 
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different n should cross. However, this effect is too small 
in magnitude to affect our results. The functional form 
we use is a convenient way to mimic the response of a 
droplet over the range of compression we investigate. 

To describe the elastic properties of a disordered 
droplet packing, we use a model which replaces the 
droplets by soft spheres which interact with their nearest 
neighbors through central-force potentials that reflect the 
behavior of the facets. However, since we found that the 
energy curve of each contact is a function of the droplet 
coordination number, the interaction energy should be 
obtained by balancing the forces at each contact, using 
the values of n for each droplet. To make the compu- 
tation tractable, we determine the average coordination 
number, n, of the whole configuration of droplets, and use 
a single effective potential for that configuration. This 
is a reasonable approximation given the rather narrow 
distribution of n for our monodisperse systems. Thus, 
we use the SE results to define a repulsive, central-force 
inter-droplet potential, 



U(d) 



2C[{^f 
0, 



If 



(d < 2R) 
(d > 2R) 



(1) 



where d = 2h is the distance between the droplet centers, 
and the factor 2 accounts for the two facets on the inter- 
acting pair. During the simulation, the exponent a(n) 
and the prefactor C{n) are estimated from cubic inter- 
polations of the values shown in Fig.|| The value of n is 
observed to increase from n c ss 6 at tp c w 0.64 to n « 10 
at ip = 0.84. We finally note that central-force potentials 
can only include compressional distortion; nevertheless 
we use them to describe shear distortions as well. 

The relaxation algorithm entails the minimization of 
the 3N-dimensional function, 
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E N = J2 U(dij), 



(2) 



i,j>i 



where dij is the distance between point-particles i and 
j. We use a conjugate-gradient (CG) method coupled 
to Brent line- minimization method [|15fl . This process 
represents a damped relaxation of a system of interact- 
ing droplets, and therefore the masses of the particles 
are not necessary. In order to find only a local mini- 
mum, our version of CG is written so that each line- 
minimization search interval is defined using the result 
of the line- minimization along the previous conjugate di- 
rection. The evaluation of the potential benefits from 
techniques borrowed from molecular dynamics, such as a 
Verlet table |L6| . The algorithm is designed to return af- 
ter a minimum number of iterations and when new values 
for the energy differ by less than 5En/En < 10~ 7 . For 
each relaxed configuration, the energy, the average coor- 
dination number, and the pressure are measured. The 
coordination number is derived from the sum of all pairs 



contributing at least a small e to the energy and val- 
ues for the pressure in each Cartesian direction (diago- 
nal elements of the stress tensor) are obtained from the 
virial |l6[ . 

We first construct a random distribution of N 
monodisperse particles of radius R in a cubic container 
with periodic boundary conditions. Typical runs start 
from a configuration prepared at the desired volume frac- 
tion. Systems that were slowly compressed from ip < <p c 
to ip > (p c , and random configurations built and relaxed 
at ip > (p c were found to have similar elastic properties. 
Thus, we usually start at ip « 0.67 > ip c in order to 
avoid the long computational time required to relax a 
configuration at p c . The system is then compressed and 
relaxed in small step increments up to ip « 0.84. From 
this value, the cubic container is sheared using isochoric 
uniaxial strains wherein, say, the z axis is stretched by a 
factor A > 1 and the perpendicular plane is compressed 
by A~^. The shear modulus is obtained from the excess 
energy density as a function of the extension ratio A [ p"7| , 

E N /V = E° N /V +\G{\ 2 + 2/A-3) (3) 

where E^ is the excess energy of the unstrained system. 
We perform several strain cycles to allow for relaxation 
and to verify reproducibility. This shearing process is 
performed along the three Cartesian directions. The sys- 
tem is then expanded in a small step increment, relaxed, 
and the shearing procedure is repeated. System sizes of 
at least N = 10 3 were used to avoid undesired relaxation 
into an fee. structure at large strains. 

The calculated ip dependence of the shear modulus 
is shown by the plus symbols in Fig. |l|, while the os- 
motic pressure is shown by the solid line, for a system 
of N — 10 3 . Interestingly, the data for n shows self- 
averaging and contains much less fluctuations that that 
for G. This suggests that there is an effective larger cor- 
relation length that determines the modulus and limits 
the averaging possible in our finite systems. Neverthe- 
less, remarkably good agreement with the experimental 
data is obtained. Both the magnitudes and the ip depen- 
dencies are correctly reproduced for both G and n. 

The simulation also provides physical insight about the 
origin of the behavior of the shear modulus of emulsions. 
There are two essential effects. The first is the positional 
relaxation of the droplets. We illustrate this graphically 
in Fig. |^: we subtract the actual motion of the droplets 
from the affinc motion caused by a strain associated with 
A = 1.0006, and use arrows at the center of each droplet 
to signify the direction and magnitude of this difference. 
The length of each arrow has been increased by 250 to 
make the motion with such a small strain observable. 
The droplet motion is clearly not afiine; moreover, the 
difference appears random in direction. The importance 
of this non-affine motion is reinforced by calculating the 
shear modulus for a strictly affine (unrelaxed) motion; it 
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FIG. 3. An isochoric uniaxial shear strain of A = 1.0006 
is applied to a TV = 10 3 system. Each arrow, magnified 
250 times, represents the droplet relaxation from an affinc, 
non-relaxed strain, tp — 0.80. 

increases by a factor of about 3 for all <p. We empha- 
size, however, that the non-affine motion of the droplets 
does not result in large scale rearrangements of their po- 
sitions. We found that only a very small fraction of 
the droplets change neighbors during shear. Thus, the 
non-affine motion results from localized relaxation of the 
droplet positions fllS]] . The second essential feature is 
the anharmonicity of the potential. We performed simi- 
lar simulations using a harmonic potential and obtained 
qualitatively different results [M] for G, which exhibited 
a significantly steeper rise near (p c , similar to what is ob- 
served in harmonic two-dimensional systems By 
contrast, the behavior of II, is not as sensitive to the 
potential. 

Finally, we note that the average coordination number 
was found to increase continuously above ip c , consistent 
with a power law increase (n — n c ) ~ (ip — ip c )^ sur- 
prisingly the same functional form as the one observed 
in computer simulations in two dimensions |P2] ]. In the 
present picture, the increase of n not only creates more 
contacts capable of storing energy, but also causes the 
energy of the existing contacts to increase because of the 
increase of the number of facets around each droplet. 
However, this has only a minor effect as confirmed by 
simulations using fixed C and a. 

These results suggest an explanation for the origin of 
the surprising elasticity of emulsions; moreover, they will 
likely have broader significance. They are immediately 
applicable to foams, since the compressibility of the gas 
in the bubbles is typically much less than the Laplace 
pressure. More generally, emulsions are an example of a 
material with strictly repulsive interactions which is nev- 
ertheless a solid when confined by an osmotic pressure. 
The elasticity of these packings is significantly different 



than that of more traditional materials, and the results 
presented here should form the basis for developing a 
more comprehensive description of these fascinating ma- 
terials. 
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